//O(N) prime sieve,
bool isprime[n];
int lp[n+1];
void sieve(int n, vi &pr) {
    pr.clear();
    fori(i,2,n) {
        if (lp[i]==0){
            pr.pb(lp[i]=i);
            isprime[i]=true;
        }
        for (int j=0;j<(int)sz(pr) && pr[j]<=lp[i] && i*pr[j]<=n;j++)
            lp[i*pr[j]]=pr[j];
    }
}

//factorization using sieve,
int lp[n+1];
void sievefac(int n, vpii &a) {
    a.clear();
    while (n>1) {
        a.pb(mp(lp[n], 1));
        if (sz(a)>1 && a[(int)sz(a)-1].first==a[(int)sz(a)-2].first) {
            a.pop_back();
            a.back().second++;
        }
        n/=lp[n];
    }
}

//Calculate phi for one number, and get factorization list
int phifac(int n, vpii &a){
    a.clear();
    int res=n;
	for (int i=2; i*i<=n; ++i)
		if (n%i==0) {
            a.pb(mp(i,0));
			for (; !(n % i); n/=i) a.back().second++;
			res-=res/i;
		}
	if (n>1) {
		res-=res/n;
		a.pb(mp(n,1));
	}
	return res;
}

//Calculate phi(i) up to N
int phi[n];
void calEulerPhi(int n){
    vi pr;
    bool check[n];
	clear0(check);
	phi[1]=1;
    fori(i,2,n) {
        if (!check[i]){
            pr.pb(i);
            phi[i]=i-1;
        }
        for (int j=0; j<(int)sz(pr)&& i*pr[j]<=n; j++){
            check[i*pr[j]]=1;
            if (i%pr[j] == 0){
                phi[i*pr[j]]=phi[i]*pr[j]; break;
            }
            else phi[i*pr[j]]=phi[i]*(pr[j]-1);
        }
    }
}

//extended euclid alogirthm
int extendgcd(int a, int b, int &x, int &y){
    if (b==0){
        x=1; y=0;
        return a;
    }
    int g=extendgcd(b, a%b, y, x);
    y-=x*(a/b);
    return g;
}

int inversemod(int a, int mod) {
    int x,y; int g=extendgcd(a,mod,x,y);
    if (g>1) return -1; return (x%mod+mod)%mod;
}

//diophatine one variable with 0<=x<mod
int diophantine1(int a, int b, int mod) {
    int g=gcd(a,mod);
    if (b%g != 0) return -1;
    a/=g; b/=g; mod/=g;
    return ((b*inversemod(a,mod))%mod);
}
//diophantine two variables with minx<=x<=maxx and miny<=y<=maxy
bool diophantinecheck (int a, int b, int c, int &x0, int &y0, int &g) {
	g=extendgcd(abs(a), abs(b), x0, y0);
	if (c%g != 0) return false;
	x0*=c/g; y0*=c/g;
	if (a < 0) x0=-x0; if (b < 0) y0=-y0;
	return true;
}

void shift_solution (int & x, int & y, int a, int b, int cnt) {
	x+=cnt*b; y-=cnt*a;
}

int diophantinefind (int a, int b, int c, int minx, int maxx, int miny, int maxy) {
	int x,y,g;
	if (!diophantinecheck(a,b,c,x,y,g)) return 0;
	a/=g; b/=g;
	int signa=(a>0)?1:-1;
	int signn=(b>0)?1:-1;
	shift_solution(x,y,a,b,(minx-x)/b);
	if (x<minx) shift_solution(x,y,a,b,signb);
	if (x>maxx) return 0;
	int lx1=x;
	shift_solution(x,y,a,b,(maxx-x)/b);
	if (x>maxx) shift_solution(x,y,a,b,-signb);
	int rx1=x;
	shift_solution(x,y,a,b,-(miny-y)/a);
	if (y<miny) shift_solution(x,y,a,b,-signa);
	if (y>maxy) return 0;
	int lx2=x;
	shift_solution(x,y,a,b,-(maxy-y)/a);
	if (y>maxy) shift_solution (x,y,a,b,signa);
	int rx2=x;
	if (lx2>rx2) swap(lx2, rx2);
	int lx=max(lx1,lx2);
	int rx=min(rx1,rx2);
	return (rx-lx)/abs(b)+1;
}

//Chinese remainder theorem
//Find x that x=remainder[i]%modu[i], all modulo must be pairwise coprime
int CRT(vi &remainder, vi &modu) {
    int productmodu=1,kq=0;
    fora(i,modu) productmodu*=modulo[i];
    fora(i,remainder){
        int x,y; extendgcd(productmodu/modu[i],modu[i],x,y);
        kq=(kq+remainder[i]*(productmodu/modu[i])*x)%productmodu;
    }
    return kq;
}

//Find x that x=remainder[i]%modu[i], no constraint
int nonCRT(vi &remainder,vi &modu){
    int ans=-1;
    fora(i,remainder) {
        if (i+1 == sz(remainder)) {
            ans=remainder[i]%modu[i]; break;
        }
        else {
            int x,y;
            int g=extend_euclid(modu[i],modu[i+1],x,y);
            if (g==0){
                ans=-1; break;
            }
            int j=remainder[i+1]-remainder[i];
            if (j%g != 0){
                ans=-1; break;
            }
            int lcm=modu[i]/g*modu[i+1];
            int d=modu[i+1]/g;
            int k=((j/g*x)%d+d)%d;
            k=(remainder[i]%lcm+k*modu[i]%lcm)%lcm;
            if (k<0) k=k+((-k)/lcm+1)*lcm;
            if (k>0) k=k-k/lcm*lcm;
            remainder[i+1]=k;
            modu[i+1]=lcm;
        }
    }
    return ans;
}

//Fast fourier transform
typedef complex<double> base;
int rev[MAXN];
//run this first as init
void fftrev(int n) {
    int logn = int(log(n)/log(2)+1);
	forn(i,n) {
		rev[i]=0;
		forn (j,logn)
			if (i & (1<<j))
				rev[i]|=1<<(logn-1-j);
	}
}

void fftmain(vector<base> &a, bool invert) {
    int n=sz(a);
	forn(i,n) if (i<rev[i]) swap(a[i],a[rev[i]]);
    for (int len=2; len<=n; len<<=1) {
		double ang=2*PI/len*(invert?-1:1);
		base wlen(cos(ang),sin(ang));
		for(int i=0; i<n; i+=len) {
			base w(1);
			for (int j=0; j<(len>>1); ++j) {
				base u=a[i+j],v=a[i+j+(len>>1)]*w;
				a[i+j]=u+v;
				a[i+j+(len>>1)]=u-v;
				w*=wlen;
			}
		}
    }
	if (invert) forn(i,n) a[i]/= n;
}

//Basic Gauss return value:
//No solution=0, One solution=1, infinity solutions = number of independent variables.
int gauss (vector < vector<double> > a, vector<double> &x) {
	int n=sz(a);
	int m=sz(a[0])-1;
	vi where (m,-1);
	for(int col=0, row=0;col<m && row<n;col++) {
		int choose=row;
		fori(i,row,n-1) if (abs(a[i][col])>abs(a[choose][col]))
            choose=i;
		if (abs(a[choose][col])<eps) continue;
		fori (i,col,m) swap(a[choose][i],a[row][i]);
		where[col]=row;
		forn(i,n) if (i != row) {
            double c=a[i][col]/a[row][col];
            fori (j,col,m) a[i][j]-=a[row][j]*c;
        }
		++row;
	}
	x.assign (m, 0);
	forn(i,m) if (where[i] != -1)
        x[i]=a[where[i]][m]/a[where[i]][i];
	forn(i,n) {
		double sum=0;
		forn(j,m) sum+=x[j]*a[i][j];
		if (abs(sum-a[i][m])>eps) return 0;
	}
    int indepen=0;
	forn(i,m) if (where[i]==-1) ++indepen;
	if (indepen>0) return 1; else return indepen;
}

// gauss modulo 2, faster version.
// return value is number of solution.
int M; //number of x
int gaussmod2 (vector < bitset<M> > a, bitset<M> &x) {
	int n=sz(a);
	int m=sz(a[0])-1;
	vi where (m,-1);
	for(int col=0, row=0;col<m && row<n;col++) {
		fori(i,row,n-1) if (a[i][col]) {
            swap(a[i],a[row]);
            break;
		}
		if (!a[row][col]) continue;
		where[col]=row;
		forn (i,n) if (i!=row && a[i][col])
            a[i]^=a[row];
		++row;
	}
	x.reset();
	forn(i,m) if (where[i] != -1) x[i]=a[where[i]][m];
	forn(i,n) {
		int sum=0;
		forn(j,m) sum=(sum+x[j]*a[i][j])&1;
		if (sum==a[i][m]) return 0;
	}
    int indepen=0;
	forn(i,m) if (where[i]==-1) ++indepen;
	return two(indepen);
}




